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Abstract 

m 

<3\ Employing a particularly suitable higher order symplectic integration 

£h algorithm, we integrate the 1-d nonlinear Schrodinger equation numerically 

| ^ for solitons moving in external potentials. In particular, we study the 

scattering off an interface separating two regions of constant potential. We 
£C) find that the soliton can break up into two solitons, eventually accompanied 

by radiation of non-solitary waves. Reflection coefficients and inelasticities 
i— ( are computed as functions of the height of the potential step and of its 

steepness. 

o 
m 

On 

^ 1 Introduction 

V 

q Recent years have seen a considerable growth in the interest for nonlinear par- 

tial differential equations with soliton solutions. In particular, the nonlinear 

q Schrodinger equation (NLSE) and its variants appear in problems drawn from 

disciplines as diverse as optics, solid state, particle and plasma physics. There, 
the NLSE describes phenomena such as modulational instability of water waves 
[1], propagation of heat pulses in anharmonic crystals, helical motion of very 
thin vortex filaments, nonlinear modulation of collisionless plasma waves [2], and 
self-trapping of light beams in optically nonlinear media [3, 4, 5]. In all these 
problems, the main interest is in the fact that the NLSE has soliton solutions. 
These are solitary waves with well defined pulse-like shapes and remarkable sta- 
bility properties [6]. 

A great deal of current interest is directed to the question how these states 
behave under the influence of external perturbations. These can be of various 
forms. We shall limit ourselves to such perturbations which can be described by 
potentials. They preserve the hamiltonian structure of the NLSE [2], but not its 
complete integrability. Other types of perturbation which also are hamiltonian 



1 



are obtained when either the coefficient in the kinetic term (the 'mass' in the 
quantum mechanical interpretation) or in the nonlinear term are made spatially 
not constant. Such inhomogeneities have indeed been studied more intensely than 
the ones we shall study below, since they are more relevant for the transmission 
of pulses through junctions in optical fibers [4, 5, 7, 8, 9]. 

More precisely, we shall consider only potentials which are constant outside 
a finite interval (we shall only consider the case of one spatial dimension). But 
we allow different values V± for x — > ±oo, mimicking thereby the effect of an 
interface between two media in which the solitons have different characteristics. 
We study initial conditions consisting of one single soliton. In general, we have 
to expect that this soliton will not just be transmitted or reflected. There might 
be also inelastic scatterings where it breaks up either into several solitons or into 
non-solitary waves, or both. 

This problem has been studied previously by several authors. While perturba- 
tive approaches were used in [3, 4, 5, 10], straightforward numerical integrations 
were made in [11]. Both approaches showed that the soliton behaves just like a 
classical particle if the force created by the potential is sufficiently weak. This 
is to be expected, but the problem what happens when the force is strong was 
left open for a simple potential ramp (the potential considered in [5] was more 
complicated). 

It is one purpose of the present paper to close this gap by means of simulations. 

Another purpose is to show the usefulness of higher order symplectic inte- 
gration algorithms. As we have already mentioned, the NLSE is a hamiltonian 
system. Thus, it is natural to apply to it integration routines which were devel- 
oped during the recent years and whose main characteristic is that they preserve 
the hamiltonian structure [12, 13, 14]. The latter is not true e.g. for standard 
methods as e.g. Runge-Kutta or predictor-corrector. Such 'symplectic' integra- 
tors (the simplest of which is the well known Verlet or 'leap frog' algorithm) have 
been applied already to the linear [15, 16, 17, 18] and nonlinear [19, 20, 21, 22] 
Schrddinger equations. 

The most popular algorithms of this type are split-operator methods. They 
depend on the hamiltonian being a sum of two terms A and B, each of which can 
be integrated explicitly. Then one uses the Baker-Campbell-Hausdorff theorem 
to approximate e t ( A + B ) t by a product of factors e takAt and e t/3kBt , where a& and (3k 
are real numbers satisfying among others aj. = (3k = 1. The error is then 
given by higher order commutators of A and B. We shall in particular apply 
a fourth order method due to McLachlan and Atela [23] which is applicable if 
one of the third order commutators vanished identically. We shall see that this 
method should be applicable to our problem, and that it is indeed numerically 
very precise, indicating that the McLachlan- Atela method is the method of choice 
for a wide class of problems. 
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2 The NLSE soliton solution 



Using appropriate units, we can write the NLSE as 
.d$(x,t) ld 2 ^{x,t) 



dt 



dx 2 



V(x,t)\ 2 V(x,t) + V(x) V(x,t), 



where V(x) is the external potential. We shall use for the latter a piecewise linear 
ansatz, with V(x) = for x < 0, V(x) = Vq > for x > xq > 0, and linearly 
rising for x between and xq, 



V(x) 




x < 

< X < Xq 
X > Xq 



(2) 



We call the negative :r-axis region /, while region II is the region x > (where 
V{x) > 0). 

We study scattering solutions where the incoming wave consists of a single 
soliton arriving from region /. The outgoing wave will then in general be a 
complicated superposition of solitons and non-solitary waves, in general moving 
both into regions / and i7. The interesting questions are how many solitons will 
leave the scattering region and with what energies, how much of the total energy 
is transmitted and reflected, and how much of it goes into non-solitary waves. 

For a constant potential Vq the soliton solutions of eq.(l) form a two-parameter 
manifold (apart from translations). Taking as parameters the velocity v and the 
amplitude a, these solutions read [24] 



cosh[a(x — vt)] 



J{vx+[(a 2 -v 2 )/2-V ]t} 



(3) 



We denote the velocity of the incoming soliton as vq. Using a suitable rescaling 
of x,t and ^, we can always choose its amplitude as ao = 1/2, without loss of 
generality. 

Among the infinitely many conserved quantities (for V(x) = const!) the 
following three are of particular interest: 
the normalization 



N = / dx 



the energy 



E 









dx 



lm 4 + v(x) 



and the momentum 



P = 



2^ 



dx 



a** 

dx 



1*1 



dx . 



dx 



(4) 



(5) 



(6) 



For the soliton given by eq.(3), N = 2a, P = vN, and E = (v 2 /2-a 2 /6)N+{V)N , 
where the average over V(x) is taken with weight oc \^\ 2 as indicated by eq.(5). 
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For a slowly varying V(x) (which implies xq/Vq ^ 1 in our case) the amplitude 
is approximately constant, and the soliton moves like a classical particle with mass 
m = 2a in an external potential mV(x) [10]. The mass of the incoming soliton 
is rriQ = 1 with our normalization. Another limit case where the soliton behaves 
like a particle is that of Vq <C Kq where Kq = Vq/2 is the kinetic energy of the 
incoming soliton. 

It is easily seen that N and E are also conserved for non-constant potential 
V, while this is not true for P. Denoting by N{, % = the normalization in 

region i, we have thus A/ iOUt + A// iOUt = Ni^ n = 1. Similarly, energy conservation 
gives £V,out + ^7/,out = E I)in = (vl - 1/12) /2. 

Conservation of N and E poses restrictions on the final state. In general, they 
do not seem to be very stringent. Assume e.g. that the final state consists of 
two solitons moving in opposite directions, (a, v) moving into region / and (6, w) 
moving into II. Then we find that 

a + b = 1/2, vl = ab + 2(av 2 + bw 2 + 2bV ) (7) 

This does not imply, in particular, a lower bound on vq since b and v can be 
arbitrarily small. Similarly, for any initial soliton we can have any number of 
outgoing solitons, provided there is at least one reflected and one transmitted 
soliton. Conservation of N and E is more stringent if no or all solitons are 
reflected. For instance, if the final state consists of a single transmitted soliton, 
then its velocity is vjj jOUt = \A>q — 2Vq. This conforms with the general statement 
that the soliton behaves like a classical particle with m = 1, and shows that there 
is no transmission if vq < \/2Vq (i.e., Kq < Vq) and xq 3> Vq. It was verified 
numerically in [11]. These authors concluded indeed that solitons impinging on 
a potential step behave like classical particles. It was mainly this claim which 
stimulated our investigation. 



3 Symplectic integration 

The NLSE is a classical hamiltonian system with Poisson bracket 

{V*(x),V(y)} = i6(x-y) (8) 

and hamiltonian H = E. This implies in particular that it can be written as 

* = = HV , (9) 

where the linear ('Liouville') operator H is defined as H ■ = { • ,H}. Split- 
operator methods can be applied by splitting H = T + V, where T and V are the 
Liouvilleans corresponding to ^ j dx\d x ^>\ 2 and J dx(— ^\^>\ 4 + Vl'I'l 2 ), 

T^= l -d 2 x ^, W = - V^f) . (10) 
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In a paper by McLachlan & Atela [23], a fourth order algorithm was introduced 
which minimizes the neglected fifth order terms in the Baker-Campbell-Hausdorff 
formula for hamiltonians for which 

{{{T,V},V},V} = 0. (11) 

This applies obviously to hamiltonians with T = \(p, M~ 1 p), V = V(x), with M 
a constant mass matrix and {qi,Pk} = $ik, since there each commutator with V 
acts as a derivative operator on any function of p. In [17] it was shown that this 
algorithm can also be applied to the linear SE where it gives better performance 
than the general fourth order algorithm [12] which does not take into account 
this special structure. 

Although the argument is less straightforward in the present case, it is not 
too hard to see that eq.(ll) holds also there [22]. Let /(|^| 2 , x) and gd^ 2 , x) be 
arbitrary functions with finite first and second derivatives. Then one finds 

{$dx\d x y\ 2 gm 2 ,x),$dyfm 2 ,y)} =ifdx(** xx *-W xx )gf, (12) 

where /' = df /d\^ 2 \, and 

rr, , * * x r , ,-. /",. ^dgd^ (x)\ 2 , x) df (\^> (x)\ 2 , x) 

UM*:** ~ ***~)<7, Jdyf} = -2. J dx\*\ 2 y[l 'J 1 ) Ml [ d J l ) . 

' ' ^ _ ' (13) 

Since the last expression is a functional of |^| 2 only, its Poisson bracket with 
J dyf vanishes identically, QED. 

The coefficients aj. and j3k for the McLachlan- Atela method are listed in [23, 
17]. Our implementation involves a spatial grid with Fourier transformation after 
each half step [17]. 

Since T and V both conserve the normalization exactly, N should be conserved 
up to round-off errors. This was checked numerically, relative errors typically were 
of order 10 -11 . Energy is not conserved exactly, and its error was fa 10~ 5 after 
an evolution time t = 300 with an integration step At = 0.005. The precise 
value depended of course on the parameters of the soliton and on xq. It was 
checked that the algorithm is indeed fourth order, and is more precise than the 
general fourth order symplectic [12] and the leap-frog (second order symplectic) 
algorithms. We also tested two other discrete Hamiltonian integration schemes 
which where examined in [26]. They both show the same qualitative behavior, 
but the discretization of the Laplace operator requires smaller time steps for the 
same spatial discretization width. All this demonstrates the advantage of the 
McLachlan-Atela algorithm. 

4 Results 

During the simulations we measured normalization Ni, energy Ei, and momentum 
Pi in each region (i = I, II) separately. The derivatives of ^> and ^* were of 
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course computed in Fourier space, as this is much more precise than taking finite 
differences in x-space. 

Since we have two conserved quantities, we can define two sets of transmission 
and reflection coefficients. We call them Tjv,i?jv and Te,Re, 

T ^ = ^T> R N = ^ = l-T N (14) 

and 

T * = lf> Re = ^ = 1~T e . (15) 

In addition we registered all local maxima of ^(rc)! 2 with ^(rc)! 2 > 1/3000. 

Since our model involves 3 free parameters (Vq,xq,vq), it is impossible to 
present results exhaustively. We did a large number of simulations with different 
parameter values, but we present only a few of them here to illustrate the variety 
of the scenarios. 

Our numerical simulations confirmed the prediction that the soliton behaves 
as a classical particle if xq/Vq 3> 1, and if Vq <C Kq. The same is true also if 
xq = and Vq = oo, i.e. if the potential acts like a hard wall. In that case, an 
exact solution of the NLSE with boundary condition \&|a;=o = and correct initial 
conditions in region I is provided by a state with two (interacting) solitons with 
opposite velocities and phases but equal amplitudes [11]. 

While the above essentially just checked the correctness of our integration 
routine, a less trivial result is that we confirmed the observations of Nogami and 
Toyama [11] for their parameter choice xq = 0, vq = 0.2, Vq ~ Kq. But we did not 
verify their claim that this is the typical behavior. Instead, the soliton typically 
breaks up and does not behave like a classical particle. 

In general, after the soliton hits the potential ramp, we found typically more 
than a single maximum of Moreover, the heights of these maxima in 

general were not constant in time, though they moved with practically constant 
velocities (see figs. 1, 3, 5, 7). Instead, they showed often very marked oscillations 
(fig. 2, 4, 6, 8) which were damped in all cases. Such damped oscillations result 
typically from superpositions of solitons with non-solitary waves [27]. We checked 
that a superposition of a soliton with a Gaussian wave packet gave essentially the 
same patterns. 

In the following we shall only show results for vq = 0.8 although, as we said, 
we had made runs also with different vq and with similar results in general. 

Figures 1 and 2 show the case where the potential is a step function (xq = 0) 
and the kinetic energy (Kq = 0.32) is larger than its height Vq = 0.3. Classi- 
cally one would expect the soliton to move into region / and to propagate there 
with a reduced speed. But our simulation shows that it breaks up into two soli- 
tons with roughly equal heights and with velocities v = —0.588 and w = 0.395. 
About half of the normalization and three thirds of the energy are transmitted 
(Tjv = 0.527, Te = 0.712). Inserting these numbers into eq.(7), we find perfect 
agreement (discrepancies are < 1%). This indicates that radiation in form of non- 
solitary waves is small in spite of the wiggles seen in fig. 2. More precisely, we 
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compared our data with ref.[27] by assuming that the transmitted wave is a single 
solitary wave immediately after leaving the interaction region. We found perfect 
agreement if we assume that this wave has exactly the same shape and width as 
the incoming soliton, but an amplitude reduced by a factor 0.728. Thus, at least 
for these parameter values, the main effect of the interaction on the transmitted 
wave is simply a reduction of amplitude. 

The situation where the potential step (Vo = 0.34, xq = 0) is higher than 
the kinetic energy Kq is plotted in the figures 3 and 4. Here one would expect 
classically that the incident soliton is completely reflected back into region II. 
But once again the behavior is quite different, the soliton splits up into two. The 
transmitted one is not as high as the reflected one (Tjv = 0.373) and therefore 
much wider, but it still carries more than half of the initial energy, Tg = 0.571. 
As we increase Vo further, the transmitted soliton rapidly shrinks. It becomes 
unobservable at Vq « 2Kq, where the soliton is practically completely reflected. 

Let us now study positive values of xq, i.e. potential ramps with finite slopes. 
Our data show unambiguously that this slope has a strong influence. If xq is 
of order 1, the soliton still breaks up as described above (figs. 5, 6), with even 
larger oscillations and even more "dirt" than for xq = 0. Flattening the potential 
ramp further but leaving its height constant, the soliton finally travels along the 
classically expected trajectory (figs. 7, 8): in the ramp region it sees a constant 
force and hence moves on a parabola; it is reflected (transmitted) for Vq > Kq 
(Vo < K ). 

This dependence on the slope of the ramp is seen very clearly when plotting 
the energy in region II as a function of time, see fig. 9. While the asymptotic 
state is reached very quickly for steep potentials, this evolution takes very long 
for gentle slopes. If xq ^> 1 (corresponding to a width of the soliton <C xq), the 
energy change is sudden when the soliton crosses the point x = 0. 

Finally, the dependence of the transmission coefficients on xq are shown in 
fig. 10. We see that they are not monotonic, with the nonmonotonicity more 
pronounced for Tjy than for Tg. This is an unexpected effect which we do not 
know how to explain. The fact that Tjv < Tg for all xq is less surprising. 

5 Summary and conclusions 

In this note we have applied an optimized fourth order symplectic integrator 
to the scattering of NLSE solitons from an external potential. The integrator 
is optimized in the sense that it takes into account that the kinetic energy is 
bilinear in It was found to be more precise than the general fourth order 
symplectic integrator. We found that solitons break in general up when hitting a 
potential threshold, in contrast to recent claims. The complexity of the outgoing 
state depends on the parameters of the potential and of the soliton, but most 
frequently the soliton breaks into two, with rather little radiation. 

The NLSE can be considered as a special case of the complex Ginzburg- 
Landau (CGL) equation $ = + a|$| 2 $ + /3V 2 $ (fJ.,a,(3 G C) with complex 
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constants. The applicability of our integrator does not depend on the phases of 
these terms, whence it should be applicable also to the CLG equation in general. 
We just have to take into account that |$| is not constant dusssring the evolution 
under the nonlinear term if Re /i, a ^ 0. In that case the integration of V involves 
solving the easy differential equation d\Q\ 2 /dt = 2(Re /i |$| 2 + Re a |$| 4 ). 

This work was partly supported by DFG within the Graduiertenkolleg "Feld- 
theoretische und numerische Methoden in der Elementarteilchen- und Statistis- 
chen Physik" , and within SFB 237. 

References 

[1] H. Hasimoto and H. Ono, J. Phys. Soc. Japan 52 (1983) 4129 

[2] G.L. Lamb, Jr., Elements of Soliton Theory (Wiley, Toronto, 1980) 

[3] Y.S. Kivshar and B.A. Malomed, Rev. Mod. Phys. 61 (1989) 763 

[4] Y.S. Kivshar, A.M. Kosevich and O.A. Chubykalo, Phys. Rev. A 41 (1990) 
1677 

[5] Y.S. Kivshar and M.L. Quiroga-Teixeiro, Phys. Rev. A 48 (1993) 4750 

[6] R. Jackiw, Rev. Mod. Phys. 49 (1977) 681 

[7] J. P. Gordon, J. Opt. Soc. Am. B 9 (1992) 91 

[8] M. Chbat et ai, J. Opt. Soc. Am. B 10 (1993) 1386 

[9] D. Anderson, M. Lisak, B. Malomed, and M. Quiroga-Teixeiro, J. Opt. Soc. 
Am. B 11 (1994) 2380 

[10] M.A. de Moura, J. Phys. A 27 (1994) 7157 

[11] Y. Nogami and F.M. Toyama, Phys. Lett. A 184 (1994) 245 

[12] H. Yoshida, Phys. Lett. A 150 (1990) 262 

[13] H. Yoshida, Celest. Mech. 56 (1993) 27 

[14] J.M. Sanz-Serna, Physica D 60 (1992) 293 

[15] A.D. Bandrauk and H. Shen, Chem. Phys. Lett. 176 (1991) 428 

[16] K. Takahashi and K. Ikeda, Kyoto preprint (1992), submitted to J. Chem. 
Phys. 

[17] H. Frauenkron and P. Grassberger, Int. J. Mod. Phys. C 5 (1994) 37 
[18] A. Rouhi and J. Wright, Comp. Phys. Commun. 85 (1995) 18 



8 



[19] D.Pathria and J.L. Morris, J. Comput. Phys. 87 (1990) 108 

[20] J.A.C. Weidemann and B.M. Herbst, SIAM J. Numer. Anal. 23 (1986) 485 

[21] A.D. Bandrauk and H. Shen, J. Phys. A 27 (1994) 7147 

[22] R.I. McLachlan, Numer. Math. 66 (1994) 465 

[23] R.I. McLachlan and P. Atela, Nonlinearity 5 (1992) 541 

[24] P.G. Drazin and R.S. Johnson, Solitons: an introduction, Cambridge Uni- 
versity Press 1989 

[25] T.R. Taha and M.J. Ablowitz, J. Computat. Phys. 55 (1984) 192, 203 
[26] M.J. Ablowitz and CM. Schober, Int. J. Mod. Phys. C 5 (1994) 397 
[27] J. Satsuma and N. Yajima, Prog. Theor. Phys. Suppl. 55 (1974) 284 



9 
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Time t 

Figure 1: Time evolution of local maxima of |^| 2 for a soliton with incident 
velocity vq = 0.8 which is scattered at a potential step with xq = and height 
Vq = 0.3 = 0.937X0. The calculation was done on a lattice with 4096 sites, 
discretization width Ax = 0.2 and integration step At = 0.005. The latter 
parameters are the same for the next figures. 
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Figure 2: Time evolution of the height of the maxima shown in fig. 1. The highest 
curve belongs to the transmitted soliton, and the second highest to the reflected 
one. The other maxima presumably are due to the superposition of non-solitary 
waves. 
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Figure 4: Same as fig. 2, but for Vq = 0.34 as in fig. 3. Now the highest curve 
belongs to the reflected soliton. 
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Figure 5: Same as fig. 1 , but for xq = 3 and Vq = 0.35 = 1.094/^0. 



T 



0.1 - 



0.01 - 



0.001 - 




100 200 300 400 500 600 700 800 

Time t 

Figure 6: Same as fig. 2, but with xq and Vq as in fig. 5. The highest curve belongs 
to the reflected soliton and the second highest to the transmitted one. The other 
maxima are side maxima presumably due to the superposition of non-solitary 
waves. 
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Figure 7: Same as fig. 1 , but for xq = 25 and Vq = 0.35 = 1.094Kq. 
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Figure 8: Same as fig. 2, but with xq and Vq as in fig. 5. 
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Figure 9: Time evolution of Ejj, i.e. the energy in region //, for different values 
of x . For all curves, v = 0.8 and V = 0.35 = 1.094A" . 
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